import itertools
import pandas as pd
import numpy as np
from docplex.mp.model import Model
import matplotlib.pyplot as plt
from haversine import haversine, Unit
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import time

start_time = time.time()  #开始时间
num_points = 111
coordinates = np.random.rand(num_points, 2)  # 生成0到1之间的随机数
# 转换为DataFrame格式，和原来的文件形式保持一致
df = pd.DataFrame(coordinates, columns=[0, 1])

t = 5
g_num = 4#产生地数
c_num = 3#收集中心数
r_num = 2#回收中心数目
fix_c = 1000 #收集中心固定费用
fix_r = 2000 #回收中心固定费用
fix_k = 1000#一级车辆的固定费用
fix_l=1000#二级车辆的固定费用
fix_z=5000#二级车辆的固定费用
T_k=100#一级车辆的运输费用
T_l=100#二级车辆的运输费用
T_z=100
T_El=1000#二级车辆的运输固废的费用
E_k=500
E_l=500
E_z=1500
lenk_c = 2
k_c=lenk_c*c_num
lenl_r=2
l_r=lenl_r*r_num
lenz_r=2
z_r=lenz_r*r_num

GM = [i for i in range (0,g_num)]#固废产生地列表
CM = [i for i in range (g_num,g_num+c_num)]#收集中心列表
RM = [i for i in range (g_num+c_num,g_num+c_num+r_num)]#回收中心列表
T = [i for i in range (t)]
T1 = [i for i in range (t-1)]

N=GM+CM+RM
CGM = CM+GM
CRM = CM+RM
GRM = GM+RM
GCRM = CM+RM+GM

#读取txt文件
x_loc = list(df.iloc[:,1])
y_loc = list(df.iloc[:,0])

#画出坐标散点图
plt.scatter(x_loc[0:g_num],y_loc[0:g_num],c='r')
plt.scatter(x_loc[g_num:g_num+c_num],y_loc[g_num:g_num+c_num],c='b',marker='d')
plt.scatter(x_loc[g_num+c_num:g_num+c_num+r_num],y_loc[g_num+c_num:g_num+c_num+r_num],c='g',marker='p')

for i in GM:
    plt.annotate('产生地/需求地''%d'%i,(x_loc[i]+0.001,y_loc[i]))
for i in CM:
    plt.annotate('一级中心''%d'%i,(x_loc[i]+0.001,y_loc[i]))
for i in RM:
    plt.annotate('二级中心''%d'%i,(x_loc[i]+0.001,y_loc[i]))
plt.show()

#各个节点之间的弧
A = [(i,j) for i in N for j in N if j!=i]

#各节点之间的距离
C = {(i, j): round(haversine(( y_loc[i],x_loc[i]), (y_loc[j],x_loc[j]), unit=Unit.KILOMETERS), 1) for i in N for j in N if j != i}

#初始化V_ijt 和 Ce_rt，作为决策变量s和x字典的键值
K=[i for i in range(k_c)]
L=[i for i in range(k_c,k_c+l_r)]
print(L)
Z=[i for i in range(k_c+l_r,k_c+l_r+z_r)]
M=K+L+Z
itm = [(i,t,m) for i in N for t in T for m in M]
mdl = Model("CDW")

#决策变量
X_i_t_m = mdl.binary_var_dict(itm,name = 'X_i_t_m')
X_c = mdl.binary_var_dict(CM,name = 'X_c')
X_r = mdl.binary_var_dict(RM,name = 'X_r')
X_k = mdl.binary_var_dict(K,name = 'X_k')
X_l = mdl.binary_var_dict(L,name = 'X_l')
X_z = mdl.binary_var_dict(Z,name = 'X_z')

#总利益
mdl.minimize(
             mdl.sum(fix_c*X_c[c] for c in CM)
             + mdl.sum(fix_k*X_k[m] for m in K)
             + mdl.sum(C[i,j]*X_i_t_m[i,t,m]*X_i_t_m[j,t+1,m]*T_k for i in CGM for j in CGM for t in T1 for m in K if j != i if j != i)
             # #承包商的利益
             + mdl.sum(fix_r * X_r[r] for r in RM)
             + mdl.sum(fix_l*X_l[m]  for m in L)
             + mdl.sum(fix_z*X_z[m] for m in Z)
             + mdl.sum(C[i, j] * X_i_t_m[i, t, m]*X_i_t_m[j,t+1,m] * T_l for i in CRM for j in CRM for t in T1 for m in L if j != i)
             + mdl.sum(C[i, j] * X_i_t_m[i, t, m] *X_i_t_m[j, t+1, m] * T_z for i in GRM for j in GRM for t in T1 for m in Z if j != i)
             # #收集中心费用
             + mdl.sum(C[i,j]*X_i_t_m[i,t,m]*X_i_t_m[j,t+1,m] * E_k for i in CGM for j in CGM for t in T1 for m in K if j != i)
             + mdl.sum(C[i, j] * X_i_t_m[i, t, m]*X_i_t_m[j,t+1,m] * E_l for i in CRM for j in CRM for t in T1 for m in L if j != i)
             + mdl.sum(C[i, j] * X_i_t_m[i, t, m] *X_i_t_m[j, t+1, m] * E_z for i in GRM for j in GRM for t in T1 for m in Z if j != i)
             # #政府要求
)
#
# 约束1确保每个一级中心和二级中心有多个车辆
for i in CM:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for m in K) <= lenk_c*X_c[i]
    )

for i in RM:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for m in L) <= lenl_r*X_r[i]
    )

for i in RM:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for m in Z) <= lenz_r*X_r[i]
    )

# 约束2确保每辆车在第1个时间步只能从一个一级中心和二级中心出发
for m in K:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for i in CM) == X_k[m]
    )

for m in L:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for i in RM) == X_l[m]
    )

for m in Z:
    mdl.add_constraint(
        mdl.sum(X_i_t_m[i, 0, m] for i in RM) == X_z[m]
    )

#约束3 确保车辆从那个厂库出回那个
for i in CM:
    for m in K:
        mdl.add_constraint(
                              mdl.sum(X_i_t_m[i,t-1,m] )
                           == mdl.sum(X_i_t_m[i,0,m])
        )

for i in RM:
    for m in L:
        mdl.add_constraint(
                              X_i_t_m[i,t-1,m]
                           == X_i_t_m[i,0,m]
        )

for i in RM:
    for m in Z:
        mdl.add_constraint(
                              mdl.sum(X_i_t_m[i,t-1,m] )
                           == mdl.sum(X_i_t_m[i,0,m])
        )

#约束4 除了二级地点，产生或开放的一级中心地点必须被访问过一次
for i in GM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,0,m] for m in K) + mdl.sum(X_i_t_m[i,0,m] for m in Z)
                       == 0
    )

for i in GM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,t-1,m]for m in K) + mdl.sum(X_i_t_m[i,t-1,m] for m in Z)
                       == 0
    )

for i in GM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,t,m] for t in T for m in K) + mdl.sum(X_i_t_m[i,t,m] for t in T for m in Z)
                       == 1
    )

for i in CM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,0,m] for m in L)
                       == 0
    )

for i in CM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,t-1,m]for m in L)
                       == 0
    )

for i in CM:
    mdl.add_constraint(
                         mdl.sum(X_i_t_m[i,t,m] for t in T for m in L)
                       == X_c[i]
    )

#约束5 车辆不能同时出现在多个地点p
for m in K:
    for t in T:
        mdl.add_constraint(
                           mdl.sum(X_i_t_m[i,t,m] for i in CGM)
                           <= X_k[m]
        )

for m in L:
    for t in T:
        mdl.add_constraint(
                           mdl.sum(X_i_t_m[i,t,m] for i in CRM)
                           <= X_l[m]
        )

for m in Z:
    for t in T:
        mdl.add_constraint(
                           mdl.sum(X_i_t_m[i,t,m] for i in GRM)
                           <= X_z[m]
        )

#约束6 确保同时步长频率正确
for m in K:
    mdl.add_constraint(
                       mdl.sum(X_i_t_m[i,t,m] for i in CGM for t in T)
                       == X_k[m]*(t+1)
    )

for m in L:
    mdl.add_constraint(
                       mdl.sum(X_i_t_m[i,t,m] for i in CRM for t in T)
                       == X_l[m]*(t+1)
    )

for m in Z:
    mdl.add_constraint(
                       mdl.sum(X_i_t_m[i,t,m] for i in GRM for t in T)
                       == X_z[m]*(t+1)
    )

solution= mdl.solve()
print(solution)

#获取决策变量的解
X_i_t_m_solution = solution.get_value_dict(X_i_t_m)
X_c_solution = solution.get_value_dict(X_c)
X_r_solution = solution.get_value_dict(X_r)

# 输出选择的收集中心和回收中心
print("收集中心选择情况:", [c for c in CM if X_c_solution[c] == 1])
print("回收中心选择情况:", [r for r in RM if X_r_solution[r] == 1])

# 获取车辆的路径
vehicle_paths = {m: [] for m in M}  # 初始化每辆车的路径
for t in T:
    for i in N:
        for m in M:
            if X_i_t_m_solution.get((i, t, m)) == 1:
                vehicle_paths[m].append(i)

# 绘制图形
plt.figure(figsize=(20, 15))

# 绘制每个节点的散点图
plt.scatter(x_loc[0:g_num], y_loc[0:g_num], c='r', label='固废产生地', s=100)
plt.scatter(x_loc[g_num:g_num+c_num], y_loc[g_num:g_num+c_num], c='b', marker='d', label='收集中心', s=100)
plt.scatter(x_loc[g_num+c_num:g_num+c_num+r_num], y_loc[g_num+c_num:g_num+c_num+r_num], c='g', marker='p', label='回收中心', s=100)

# 标注每个点
for i in GM:
    plt.annotate(f'产生地 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')
for i in CM:
    plt.annotate(f'一级中心 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')
for i in RM:
    plt.annotate(f'二级中心 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')

# 绘制每辆车的路径
for m in M:
    path_x = [x_loc[i] for i in vehicle_paths[m]]
    path_y = [y_loc[i] for i in vehicle_paths[m]]
    plt.plot(path_x, path_y, marker='o', label=f'车辆 {m}')

# 绘制图形并添加箭头表示时间顺序
plt.figure(figsize=(5, 4))

# 绘制每个节点的散点图
plt.scatter(x_loc[0:g_num], y_loc[0:g_num], c='r', label='固废产生地', s=100)
plt.scatter(x_loc[g_num:g_num + c_num], y_loc[g_num:g_num + c_num], c='b', marker='d', label='收集中心', s=100)
plt.scatter(x_loc[g_num + c_num:g_num + c_num + r_num], y_loc[g_num + c_num:g_num + c_num + r_num], c='g', marker='p',
            label='回收中心', s=100)

# 标注每个点
for i in GM:
    plt.annotate(f'产生地 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')
for i in CM:
    plt.annotate(f'一级中心 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')
for i in RM:
    plt.annotate(f'二级中心 {i}', (x_loc[i], y_loc[i]), fontsize=10, ha='right')

# 绘制每辆车的路径，并添加箭头指示时间顺序
for m in M:
    path_x = [x_loc[i] for i in vehicle_paths[m]]
    path_y = [y_loc[i] for i in vehicle_paths[m]]

    # 绘制路径
    # plt.plot(path_x, path_y, marker='o', label=f'车辆 {m}')

    # 添加箭头（从每个点到下一个点，表示车辆的行驶方向）
    for i in range(len(path_x) - 1):
        plt.arrow(path_x[i], path_y[i], path_x[i + 1] - path_x[i], path_y[i + 1] - path_y[i],
                  head_width=0.005, head_length=0.01, fc='black', ec='black')

# 设置图形属性
plt.title("优化的车辆路径（带时间顺序）")
plt.xlabel("经度")
plt.ylabel("纬度")
# plt.legend()
plt.grid(True)
plt.show()